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Application  of  Chole sky -Like  Matrix  Decomposition  Methods 
to  the  Evaluation  of  Atomic  Orbital  Integrals 
and  Integral  Derivatives 

Douglas  W.  O'Neal^  and  Jack  Simons 
Department  of  Chemistry 
University  of  Utah 
Salt  Lake  City,  UT  84112 

Abstract 

When  viewed  as  a  square  two -indexed  matrix,  the  array  of  atomic 
orbital  based,  two-electron  integrals  (ij|kl)  is  a  positive 
semidefinite  array.  Beebe  and  Linderberg  showed,  in  1977,  that  actual 
or  near  linear  dependencies  often  exist  within  the  types  of  atomic 
orbital  basis  sets  employed  in  conventional  quantum  chemical 
calculations.  In  fact,  large  (i.e.,  higher  quality)  bases  were  shown 
to  be  substantially  more  redundant  than  smaller  or  more  spatially 
separated  bases.  In  situations  where  there  exists  significant  basis 
near  redundancy,  the  rank  (r)  of  the  (ij|kl)  -  V  matrix  of 
integrals  will  be  significantly  smaller  than  the  matrix  dimension  M. 
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The  idea  of  using  the  Cholesky  procedure  for  generating 
two-electron  integrals  was  put  forth  by  Beebe  and  Linderberg  in  19771. 
Linear  dependencies  caused  by  either  basis  set  redundance  or  symmetry 
give  the  matrix  (ij|kl)  a  true  rank  (r)  smaller  than  the  apparent  rank 
(i.e.,  its  dimension  M) .  In  applying  the  method  to  integral 
derivatives  as  described  in  Section  la,  additional  dependencies  can 
arise.  The  geometrical  derivative  of  an  integral  involving 

Gaussian- type  orbitals  (GTO's)  can  be  expressed  as  a  difference 
between  integrals  over  GTO's  with  angular  momentum  quantum  numbers  one 
higher  and  one  lower  than  that  of  the  differentiated  GTO.  Thus,  if  a 
very  flexible  basis  set  is  used,  one  might  reasonably  expect  that  many 
of  the  integral  derivatives  (i.e.,  those  involving  GTO's  with  smaller 
angular  momentum)  could  be  reexpressed  as  linear  combinations  of 
two-electron  integrals  themselves,  in  which  case  redundancy  would 
occur  in  the  combined  list  of  integrals  and  integral  derivatives. 

In  the  next  two  subsections  we  describe  how  we  implement  the 
decomposition  of  the  integral  and  integral  derivative  array  and  how  we 
have  computed  those  subsets  of  integrals  and  integral  derivatives  that 
the  Cholesky  algorithm  requires.  In  Section  II  we  present  and 

analyze  the  results  of  our  calculations  on  C  and  on  benzene. 

Ia .  Procedure 

Consider  the  atomic  orbital  basis  {^(r)},  (k-l,...,N)  where  N  is 
the  number  of  basis  functions.  The  two-electron  integrals  form  a 


matrix  V  defined  as 
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V  -  V  -  (ij I kl) 

I  J  i  j  ,  k  1  J  1 

-  IT  **(l)rf.(l)  ~~  **(2)^(2)  dr  dr  (1) 

To  simplify  notation,  we  use  a  capital  index  to  denote  the  two  indices 
occurring  in  the  orbital  products  I-ij ,  J-kl.  The  range  of  I  and  J 
is  then  from  1  to  M  where  M  -  N(N+l)/2. 

The  matrix  of  two-electron  integrals  is  clearly  positive 
semi-definite  as  1/r^  is  a  positive  definite  weighting  factor  with 
the  volume  elements  (dr^dr^ .  It  can  be  decomposed  via  the  Cholesky 
algorithm  to  form  a  lower  triangular  matrix  i  such  that 

¥  -  LLt  (2) 

As  long  as  V  is  positive  definite,  L  always  exists,  and  it  has  been 
shown  that  the  Cholesky  algorithm  is  numerically  stable.2 

The  algorithm  utilized  to  decompose  V  into  can  be  detailed  as 
follows.  For  J-1,2,...,M,  we  compute 


[J  “  1  -1 

v- l  «• , 

k«l  J 

-  [  vi.J  TTaJ  r-  ■  for  W+1 . M  (3-2) 

L  k-l  J  J,J 

If  carried  through  to  completion,  in  which  case  L  would  have  the  same 
dimension  as  V  itself,  this  process  would  require  the  evaluation  of  M 
square  roots,  M(M-l) (M-4)/6  multiplications  or  divisions,  and 
M(M-l) (M+2)/6  additions  or  subtractions.  It  would  also  necessitate 
the  calculation  of  each  of  the  elements  of  the  original  V  array,  and 
would  therefore  not  produce  any  savings.  In  fact,  since  M-N(N+l)/2, 
this  process,  if  carried  to  completion,  leads  to  on  the  order  of  N6/12 
arithmetic  steps,  in  which  case  the  decomposition  of  V  into  I*  would  be 
more  time  consuming  than  the  evaluation  of  the  full  V  matrix  by 


1/2 


(3.1) 


and,  for  each  J,  we  evaluate 
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conventional  means . 

The  key  to  implementing  the  above  algorithm  in  a  manner  which  can 
overcome  this  N8  dependence  lies  in  processing  the  matrix  elements  in 
a  way  which  allows  the  (potential)  reduced  rank  of  V  to  be  exploited. 
We  emphasize  that  it  is  the  physical  content  (i.e.,  near  redundancy 
and  saturation)  of  the  basis,  which  gives  rise  to  the  rank  reduction 
and  allows  the  process  we  describe  here  to  succeed,  rather  than  the 
process  itself.  To  perform  the  Cholesky  decomposition  in  an  order 
which  exploits  rank  reduction,  we  proceed  as  follows: 

1)  We  first  calculate  all  of  the  diagonal  elements  V  using 

conventional  integral  evaluation  methods,  and  we  sort 

these  diagonal  V  into  a  non- increasing  order,  keeping  a 
J,  J 

record  of  the  original  order. 

1/2 

2)  For  J  -  1,2, .. .M,  we  set  L  -  V  ,  and  we  then  (for  each  J) 

J  ,  J  J ,  J 

3)  calculate  the  column  of  integrals  V  ,  for  I  -  J+1,...,M 

X ,  J 

using  conventional  means. 

4)  We  then  (for  each  J)  set  L  -  [v  -  7  L  L  — , 

i.j  I  i.J  L  i , k  j  , icj  L 

L  k- 1  J  J , J 

for  I  -  J+l, . . . ,M. 


5)  Finally,  we  modify  the  diagonal  elements  according  to 


V  -  V  -  L  , 

i, i  x.i  I.J 


I  -  J+l, 


.M 


This  is  equivalent  to  the  procedure  given  in  equations  3 . 1  and  3.2. 
However,  at  any  point  in  the  calculation  (say  the  Jth  step),  V  is 
the  largest  remaining  element  in  the  matrix  V.  If  the  matrix  elements 
are  known  or  needed  only  to  a  given  numerical  accuracy  5 ,  then  when 
V  3  5  the  calculation  can  be  stopped  with  no  numerical  loss.  If 
the  decomposition  is  stopped  at  point  J,  then  in  the  reformation  of 
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the  V  matrix  as  LLT,  all  elements  formed  in  steps  1  to  J  are  exact  (to 
machine  precision) ,  and  elements  formed  in  steps  J+l  to  M  are  in  error 
by  less  than  S.  Clearly  the  key  question  is  whether,  for  reasonable 

”9 

values  of  S  (e.g.,  5  ~  10  ),  the  decomposition  process  can  be  stopped 
(at  J-r)  early  enough  to  make  this  decomposition  less  time  consuming 
than  the  evaluation  of  the  full  matrix  by  conventional  means.  The 
process  itself  will  not  succeed  unless,  for  a  given  cutoff  5,  the  rank 
r  of  V  grows  slower  than  the  dimension  M  of  V;  that  is,  the  (near) 
redundancy  in  the  basis  function  space  is  what  can  make  the  process 
succeed. 

In  expanding  the  method  to  include  first  integral  derivatives, 
the  matrix  V  is  redefined,  but  the  machinery  of  the  decomposition 
given  above  is  unchanged.  When  dealing  with  the  integrals,  a  row  (or 
column)  of  V  is  labeled  by  the  functions  {f  )  -  {0*(1)0  (1),  i  -  1,N; 

i  J  i  J 

j  -  1,  i).  To  include  the  integral  derivative  case,  this  set  of 

functions  is  expanded  to  be  (f  )  -  {0*(1)0  (1)  ,^-(0*(l)0  (1) ) )  ,  where 

ij  i  j  QX  i  j 

x  labels  the  atomic  displacements  whose  integral  derivatives  are  to  be 
included.  If  the  derivative  function  is  zero,  the  corresponding  f 
is  deleted  from  the  set. 

Ib.  Integral  Evaluation  Methods 

As  was  stated  before,  the  Cholesky  decomposition  takes  advantage 
of  (near)  redundancies  in  the  matrix  of  two-electron  integrals.  The 
manner  of  formation  of  the  integrals  or  even  the  type  of  integrals 
does  not  enter  directly  into  the  procedure.  Gaussian- type  orbitals 
were  chosen  in  this  implementation.  For  a  GTO,  the  0's  in  eq  1  takes 
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the  form 

2 

*(l,m,n,a,R)  -  (x-X)1(y-Y)"(z-Z)n  e'aR  . 

The  use  of  GTO's  is  widespread  in  quantum  chemical  calculations 
because  of  the  relative  ease  with  which  electron  repulsion  integrals 
can  be  calculated  compared  with  other  basis  set  types  such  as  Slater- 
type  orbitals. 

The  Rys  polynomial  method  of  two  electron  integral  and  integral 

3 

derivative  evaluation  is  used.  This  method  was  chosen  because  of  the 

ease  and  efficiency  with  which  integrals  over  high  angular  momentum 

functions  can  be  calculated.  In  this  method,  a  Gaussian  transform  of 

the  Coulomb  operator  in  eq  1  is  performed 

r  “  Vw  J  «*P<-ur  )du.  (4) 

12  0 

This  allows  a  separation  of  variables  for  the  three  cartesian 
coordinates  so  that  the  two  electron  integral  can  be  rewritten  as 

(ij|kl)  -  J”  i;<u)r  (u)I'(u)  du.  (5) 

I'  ,  I'  ,  and  I'  are  two-dimensional  integrals  over  the  respective 
cartesian  coordinates  of  the  two  electrons.  With  a  change  of 
variables 


.2  2  .  ,  2. 
t  -  U  /  (p+U  ) 


I' 

X 


Ix  -  1//  1-t2 
p  -  (a+/3)  (y+6) /(a+ft+y+S) 
the  two -electron  integral  takes  the  form 

(ij|kl)  -  2 Jl  f1  I  (t)I  (t)I  (t)  dt 

v  /f  j  x  y  z 
0 

-  f  pl<t) 


e-/,PQ2t2(jt 


(6.1) 

(6.2) 

(6.3) 


(7) 


fch  2 

where  PL(t)  is  an  L  degree  polynomial  in  t  ,  L-l^+1^+1^+1^,  and  P 
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and  Q  are  functions  of  the  Gaussian  exponents  and  positions. 

When  one  differentiates  the  expression  in  eq  7  with  respect  to 
position  on  one  nucleus,  the  derivative  passes  through  the  integral. 

th  2 

The  integrand  is  then  an  L+l  polynomial  in  t  and  may  be  evaluated 
using  the  Rys  quadrature. 


ax 


3(lJ|kl)  r-  SI  <t) 

-  -  2/i  I  Tic—  yc>  i,<t> 3t 

0  a 

-  f1  P'  (t)  e 
J  L+l 


2  2 
-pPQ  tZ 


at 


n'ai  *  ,*1 

-  [  1  —  I_(t J  IJtJ  V_  ,  n'  >  ^ 


a-  i  ax 


y  a  x  a 


(8) 


The  recursion  relationship  given  in  reference  3c  has  been 


dlx( t  ) 

differentiated  to  give  the  quantity  — - -  directly. 

ax 

a 

The  integral  derivative  could  have  been  expressed  as  a  difference 
of  integrals  over  higher  and  lower  angular  momentum  functions,  but,  as 

4 

Pulay  has  pointed  out  ,  the  relative  efficiency  of  the  two  methods  is 
questionable.  The  current  implementation  was  judged  to  be  more  easily 
computer  coded  and  was  chosen  for  that  reason. 


II.  Results 

Ila.  The  Test  Calculations 

As  the  major  test  of  this  procedure,  a  series  of  calculations  was 
done  on  diatomic  carbon  with  twenty- two  separate,  fully  uncontracted, 
even-tempered  basis  sets  varying  in  size  from  ten  (2slp)  to 
seventy- two  (18s6p)  functions.  This  large  range  of  basis  sizes  was 
studied  to  explore  the  effects  of  increasing  basis  near  redundancy. 


In  each  eest  calculation,  the  Cholesky  decomposition  was  allowed  to 
proceed  to  tolerance  ( 6 )  of  zero.  As  the  process  advanced, 
intermediate  rank  reduction  information  (i.e.  r  as  a  function  of  S ) 
and  CPU  timing  data  were  saved.  This  was  done  first  on  the  list  of 
two-electron  integrals  and  then  for  the  list  of  integrals  and  integral 
derivatives  (with  respect  to  the  C-C  bond  length  coordinate).  Tables 
1-4  summarize  the  rank  reduction  data  realized  in  these  calculations 
at  S  tolerance  values  of  10  3 ,  10  8 ,  10  8 ,  10  12 ,  and  10  15  (in  most 
quantum  chemistry  codes,  integrals  and  integral  derivatives  are 
computed  to  precision  of  approximately  10  11  au) .  Graphical 
representations  of  this  data  are  shown  in  figures  1-4. 

All  calculations  shown  here  were  done  on  our  FPS-164  array 
processor  but  did  not  make  use  of  its  matrix  accelerator  (MAX)  board. 
For  the  purposes  of  testing  the  Cholesky  algorithm,  all  integrals 
which  were  not  zero  by  symmetry  were  explicitly  calculated  although 
this  should  not  be  necessary  in  production  runs. 

lib.  Analysis  of  Findings 

As  the  data  of  Tables  1  and  2  show,  the  rank  r  of  the  matrices  is 
much  smaller  than  the  dimension  M  even  for  tight  cutoff  tolerances  and 
increases  slowly  as  the  dimension  increases.  As  M  continues  to 
increase,  the  rank  seems  to  approach  an  asymptotic  value  (see  figs  3 
and  4)  after  which  addition  of  basis  functions  does  not  increase  the 
rank.  At  this  point,  the  basis  set  is  complete  to  within  the  cutoff 
tolerance.  To  analyze  the  potential  CPU  and/or  storage  advantages  of 
the  Cholesky  process  applied  to  the  integrals  and  integral 
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derivatives,  we  wanted  to  fit  the  "data"  obtained  in  our  calculations 
to  a  reasonable  functional  form. 

In  choosing  a  functional  form  to  fit  the  M  dependence  of  r,  we 
required  the  function  have  an  asymptotic  limit  as  the  above  analysis 
indicates  and  a  slope  of  approximately  one  for  small  values  of  M 
because  r  ->  M  in  this  limit.  We  therefore  model  the  variation  of  r 
with  M  as  atanh(/?M7) .  For  this  function,  the  asymptotic  (as  M  -»  «) 
limit  is  a  and,  if  the  optimal  value  of  &  is  approximately  1/a,  the 
slope  is  approximately  one  for  small  M. 

Results  of  least  squares  fits  of  the  above  functional  form  to  our 

r(M)  "data"  for  the  two-electron  integrals  for  diatomic  carbon  are 

given  in  table  5  and  for  two -electron  integrals  and  derivatives  in 

table  6  for  the  five  cutoff  tolerances  examined.  Figures  5  and  6  show 

the  resulting  optimal  leasts  squares  curves  for  the  two-electron 

integrals  and  for  the  integrals  and  derivatives  respectively.  In  the 

two-electron  integral  case,  the  a  values  range  from  186  to  1012  with  P 

values  between  0.01140  and  0.00066  and  7  vaues  between  1.063  and 

1.761.  For  the  two-electron  integrals  and  derivatives,  a  ranges  from 

326  to  1324,  p  from  0.01304  to  0.00110,  and  7  from  1.028  to  1.617. 

One  notes  here  that  the  values  of  p  are  approximately  1/a  for  each 

case,  satisfying  the  initial  slope  requirement  (i.e.,  the  fit  to  the 

actual  "data"  is  quite  good  as  K  -t  0) .  The  result  that  a  values  are 

less  than  2500  is  of  crucial  importance  as  we  now  illustrate 

2 

(corresponding  to  basis  set  sizes  near  50  because  H  «  N  ). 


lie.  Potential  Numerical  Advantage  of  the  Method 
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Analysis  of  the  algorithm  described  in  Section  la.  shows  that 
the  Cholesky  decomposition  process  requires:  (1)  arithmetic 

operations  (additions  or  multiplications)  whose  numbers  scale  as 
(3M-2r+l) (r-l)r/3  as  M  (and  r)  increases  and  (2)  integral  or  integral 
derivative  evaluations  whose  number  varies  as  Mr  -  r(r-l)/2.  In 
contrast,  the  conventional  treatment  requires  M(M+l)/2  integrals  to  be 
evaluated.  If  the  CPU  time  required  for  a  floating  point  addition  or 
multiplication  ij  denoted  A  and  the  (average)  time  required  to 
evaluate  a  single  two-electron  integral  or  derivative  is  denoted  B, 
the  ratio  of  the  Cholesky  to  conventional  computer  time  requirements 
should  vary  as 

r/T  A(3M-2r+l) (r-l)r/3  +  B(Mr-r(r-l)/2) 

'  "  BM(M+l)/2 

This  ratio  must  be  less  than  one  for  the  Cholesky  process  to  be 
computationally  advantageous.  On  our  FPS-164  array  processor  using 
non-vectorized  code,  a  floating  point  multiplication  takes  0.540  nsec 
(with  an  addition  being  a  factor  of  three  less)  and  a  two-electron 
integral  or  integral  derivative  evaluation  requires  503  psec  using  the 
highly  optimized  ARGOS  integral  evaluation  program6.  Thus  for  the 
following  analysis,  we  set  A-0.540  and  B-503. 

Figures  7  and  8  show  C/T  for  the  diatomic  carbon  data  for 
integrals  and  integral  plus  derivatives  respectively.  For  the 

integrals,  the  ratio  C/T  is  less  than  one  for  all  but  the  tightest 
cutoff  tolerance  and  is  a  decreasing  function  of  M  as  M  increases. 
For  the  integrals  and  derivatives,  the  result  is  more  pronounced  with 
C/T  being  less  than  0.4  for  all  values  of  M  and  decreasing  to  below 
0.1  for  M  >  1000.  Thus,  the  Cholesky  decomposition  process  shows  much 
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promise  for  decreasing  computation  times  for  integral  and 

(especially)  integral  and  derivative  evaluation. 

Disk  storage  requirements  may  also  be  greatly  reduced  by  the 
Cholesky  process.  The  conventional  algorithms  store  on  the  order  of 
M2  integrals.  The  Cholesky  algorithm  described  in  Section  la 
stores  M(r+1)  -  r(r+l)/2  integrals.  The  ratio  of  the  Cholesky  to 
conventional  storage  requirements  is  less  than  unity  for  all  M  and, 
for  large  M  (as  r  approaches  its  asymptotic  constant  a)  decreases  as 
a/M.  Saving  substantially  on  storage  also  leads  to  increased  program 
performance  through  reduced  I/O.  With  large  main  memories 
characteristic  of  modern  computers,  it  is  possible  that  all  of  the 
requisite  integrals  can  be  kept  in  the  high-speed  memory  of  the 
machine . 

In  addition  to  savings  in  integral  evaluation  time  and  storage, 
the  Cholesky  process,  if  advantageous,  can  give  major  savings  when 
atomic  orbital  based,  two-electron  integrals  (ij|kl)  must  be 
transformed  to  the  molecular  orbital  (MO)  basis.  Because  each 


integral  (those  given  exactly  and  those  approximately)  is  expressed  in 

r 

the  Cholesky  procedure  as  V  -  V  L  L  ,  the  transformation  of  the 

I .  J  “  I  ,  k  k ,  J 

k-1 


I  and  J  indices  to  the  mo  basis  set  can  be  realized  by  transforming 


the  first  (I)  index  of  the  {1^  array  to  the  mo  basis: 

r  L  CC  -L  -L 

ij.k  ia  jb  ab,k  A,k 

where  the  {C  )  are  the  LCAO-MO  expansion  coefficients. 


(9) 

This 


two- index  transformation  would  be  carried  out  for  all  k-l,...,r,  but 


the  k  index  itself  need  not  be  transformed.  Thus,  the  requisite 
transformation  process  involves  rN(N+l)  steps  to  obtain  { } ;  the 
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conventional  transformation  of  V  to  requires  4N5  steps. 

lid.  Other  Tests 

As  a  check  on  the  generality  of  our  findings  on  using  the 
above  bases ,  two  other  sets  of  calculations  were  performed  with 
substantially  different  basis  sets.  A  Dunning  10s6p  primitive  set7 
and  an  uncontracted  6-311G  set  (lls5p  primitives)  were  used  for  both 
integral  and  integral  plus  integral  derivative  calculations.  Table  7 
summarizes  the  results  for  these  calculations.  Comparing  these  results 
with  those  using  similarly  sized  even-tempered  basis  sets  shows  that 
the  results  are  very  similar.  We  therefore  feel  it  is  likely  that 
integral  and  derivative  evaluations  within  any  reasonably  well 
optimized  moderate  to  large  size  basis  set  will  benefit  from  the 
Cholesky  decomposition. 

An  additional  calculation  was  done  on  the  benzene  molecule  to 
judge  the  performance  of  the  Cholesky  decomposition  on  a  polyatomic 
system  for  which  many  integrals  vanish  due  to  large  spatial  separation 
between  atomic  centers  and  for  which  basis  set  saturation  was  not 
anticipated.  The  36  function  ST0-3G  basis  of  Hehre,  Stewart,  and 

9 

Pople  was  used  in  this  calculation  and  the  integral  derivatives  were 
taken  with  respect  to  one  C-C  bond  distance.  Although  the  rank 
reduction  is  still  significant,  the  results  (Table  8)  show  a  marked 
decline  in  the  advantage  of  the  Cholesky  process.  This  is  primarily 
due  to  the  smaller  basis  set  size  and  large  spatial  separations  among 
the  orbitals  in  the  molecule.  With  this  much  less  flexible  basis, 
there  will  be  fewer  near  dependencies  in  the  basis,  and  the  Cholesky 
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process  is  not  advantageous . 

Ill .  Summary 

The  Cholesky  decomposition  method  has  been  shown  to  be  effective 
in  reducing  the  number  of  integrals  and  integral  derivatives  that  need 
to  be  calculated,  stored,  and  transformed  to  the  molecular  orbital 
basis  for  atomic  orbital  basis  sets  which  contain  substantial  near 
redundancy.  The  saturation  of  the  atomic  orbital  basis,  not  the 
Cholesky  process  itself,  causes  this  procedure  to  succeed.  As 
analytical  energy  and  wavefunction  derivative  methods  are  applied  to 
larger  molecular  systems,  the  number  of  requisite  integrals  and 
integral  derivatives  grow  very  rapidly.  The  incorporation  of  our 
method  into  existing  analytical  derivative  program  suites  should 
therefore  allow  CPU  and  disk  storage  requirements  to  grow  at  a  much 
slower  pace  and  thereby  allow  calculations  on  much  larger  systems  than 
is  thought  currently  practical. 
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Table  1  Rank  reduction  at  cutoff  tolerances  of  10  ,  10  ,  10  , 


10  12 ,  and  10  15  for  C2  two- electron  integrals. 


basis  set 

#  basis 
functions 

full 

rank 

rank  at  rank  at  rank  at  rank  at 
10**- 3  10**- 6  10**- 9  10**- 12 

rank  at 
10**- 15 

2slp 

10 

55 

25 

31 

32 

32 

32 

3slp 

12 

78 

30 

45 

46 

46 

46 

4s  lp 

14 

105 

34 

60 

63 

64 

64 

5slp 

16 

136 

37 

72 

80 

81 

82 

5s2p 

22 

253 

60 

106 

134 

166 

145 

6s2p 

24 

300 

58 

123 

165 

174 

179 

7s2p 

26 

351 

64 

133 

185 

198 

208 

8s2p 

28 

406 

66 

138 

196 

214 

232 

8s3p 

34 

595 

87 

170 

257 

308 

328 

9s3p 

36 

666 

91 

182 

278 

334 

360 

10s3p 

38 

741 

91 

187 

288 

353 

396 

lls3p 

40 

820 

98 

194 

295 

385 

424 

lls4p 

46 

1081 

111 

230 

356 

472 

541 

12s4p 

48 

1176 

116 

234 

360 

481 

558 

13s4p 

50 

1275 

115 

240 

364 

491 

564 

14s4p 

52 

1378 

114 

244 

375 

505 

603 

14s5p 

58 

1711 

130 

276 

422 

580 

684 

15s5p 

60 

1830 

131 

275 

423 

585 

732 

16s5p 

62 

1953 

134 

280 

432 

598 

739 

17s5p 

64 

2080 

134 

285 

448 

603 

766 

17s6p 

70 

2485 

148 

319 

481 

670 

850 

18s6p 

72 

2628 

149 

324 

495 

674 

863 
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Table  2  Rank  reduction  at  cutoff  tolerances  of  10  3 ,  10  6 ,  10  9 , 
10  12 ,  and  10*15  for  C  two -electron  integrals  and 
derivatives . 


basis  set 

#  basis 
functions 

full 

rank* 

rank  at 
10**- 3 

rank  at 
10** -6 

rank  at 
10**- 9 

rank  at 
10**-12 

rank  at 
10**- 15 

2slp 

10 

110 

47 

54 

55 

55 

55 

3slp 

12 

156 

58 

77 

79 

79 

79 

4s  lp 

14 

210 

67 

99 

103 

104 

104 

5s  lp 

16 

272 

69 

112 

122 

124 

125, 

5s2p 

22 

506 

104 

159 

188 

222 

203 

6s2p 

24 

600 

103 

181 

227 

237 

243 

7s2p 

26 

702 

111 

195 

251 

265 

276 

8s2p 

28 

812 

114 

203 

264 

284 

303 

8s3p 

34 

1190 

156 

258 

356 

408 

430 

9s3p 

36 

1332 

156 

269 

373 

433 

461 

10s3p 

38 

1482 

157 

280 

390 

459 

503 

lls3p 

40 

1640 

169 

287 

399 

492 

534 

lls4p 

46 

2162 

202 

361 

500 

619 

690 

12s4p 

48 

2352 

203 

358 

496 

624 

702 

13s4p 

50 

2550 

201 

355 

493 

627 

701 

14s4p 

52 

2756 

205 

363 

511 

645 

747 

14s5p 

58 

3422 

223 

407 

574 

738 

850 

15s5p 

60 

3660 

235 

415 

585 

755 

908 

16s5p 

62 

3906 

236 

418 

593 

767 

913 

17s5p 

64 

4160 

226 

420 

610 

770 

937 

17s6p 

70 

4970 

255 

477 

669 

867 

1058 

18s6p 

72 

5256 

258 

482 

685 

872 

1070 

"includes  orbital  products  and  derivatives 
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Table  3  Fraction  of  apparent  rank  at  cutoff  tolerances  of  10  , 

10  ®,  10~9,  10  12,  and  10  15  for  C 2  two-electron  integrals. 


basis  set 

#  basis 
functions 

full 

rank 

frac.of 
rank  at 
10**- 3 

frac . of 
rank  at 
10** -6 

frac.of 
rank  at 
10**-9 

frac.of 
rank  at 
10**- 12 

frac.of 
rank  at 
10**- 15 

2slp 

10 

55 

0.45454 

0.56363 

0.58181 

0.58181 

0.58181 

3slp 

12 

78 

0.38461 

0.57692 

0.58974 

0.58974 

0.58974 

4slp 

14 

105 

0.32381 

0.57142 

0 . 60000 

0.60952 

0.60952 

Sslp 

16 

136 

0.27205 

0.52941 

0.58823 

0.59558 

0.60294 

5s2p 

22 

253 

0.23715 

0.41897 

0.52964 

0.65612 

0.57312 

6s2p 

24 

300 

0.19333 

0.41000 

0.55000 

0 . 58000 

0.59666 

7s2p 

26 

351 

0.18233 

0.37891 

0.52706 

0.56410 

0.59259 

8s2p 

28 

406 

0.16256 

0.33990 

0.48275 

0.52709 

0.57142 

8s3p 

34 

595 

0.14621 

0.28571 

0.43193 

0.51764 

0.55126 

9s3p 

36 

666 

0.13663 

0.27327 

0.41741 

0.50150 

0.54054 

10s  3p 

38 

741 

0.12280 

0.25236 

0.38866 

0.47638 

0.53441 

lls3p 

40 

820 

0.11951 

0.23658 

0.35975 

0.46951 

0.51707 

lls4p 

46 

1081 

0.10268 

0.21276 

0.32932 

0.43663 

0.50046 

12s4p 

48 

1176 

0.09863 

0.19898 

0.30612 

0.40901 

0.47449 

13s4p 

50 

1275 

0.09019 

0.18823 

0.28549 

0.38509 

0.44235 

14s4p 

52 

1378 

0.08272 

0.17706 

0.27213 

0.36647 

0.43759 

14s5p 

58 

1711 

0.07597 

0.16130 

0.24663 

0.33898 

0.39976 

15s5p 

60 

1830 

0.07158 

0.15027 

0.23114 

0.31967 

0.40000 

16s5p 

62 

1953 

0.06861 

0.14336 

0.22119 

0.30619 

0.37839 

17s5p 

64 

2080 

0.06442 

0.13701 

0.21538 

0.28990 

0.36826 

17s6p 

70 

2485 

0.05955 

0.12837 

0.19356 

0.26961 

0.34205 

18s6p 

72 

2628 

0.05669 

0.12328 

0.18835 

0.25646 

0.32838 
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Table  4  Fraction  of  apparent  rank  at  cutoff  tolerances  of  10  , 

10  8 ,  10  9 ,  10  12 ,  and  10  15  for  C  two -electron  integrals 
and  derivatives. 


basis 

set 

#  basis 
functions 

full 

rank3 

f rac . of 
rank  at 
10**- 3 

f rac .of 
rank  at 
10**- 6 

f rac. of 
rank  at 
10**- 9 

f rac. of 
rank  at 
10**- 12 

f rac . of 
rank  at 
10**- 15 

2slp 

10 

110 

0.42727 

0.49091 

0.50000 

0 . 50000 

0.50000 

3slp 

12 

156 

0.37179 

0.49359 

0.50641 

0.50641 

0.50641 

4s  lp 

14 

210 

0.31905 

0.47142 

0.49048 

0.49524 

0.49524 

5slp 

16 

272 

0.25368 

0.41176 

0.44853 

0.45588 

0.45956 

5s2p 

22 

506 

0.20553 

0.31423 

0.37154 

0.43874 

0.40119 

6s2p 

24 

600 

0.17167 

0.30167 

0.37833 

0.39500 

0.40500 

7s2p 

26 

702 

0.15812 

0.27778 

0.35755 

0.37749 

0.39316 

8s2p 

28 

812 

0.14039 

0.25000 

0.32512 

0.34975 

0.37315 

8s3p 

34 

1190 

0.13109 

0.21681 

0.29916 

0.34286 

0.36134 

9s3p 

36 

1332 

0.11712 

0.20195 

0.28003 

0.32508 

0.34610 

10s  3p 

38 

1482 

0.10594 

0.18893 

0.26316 

0.30972 

0.33941 

lls3p 

40 

1640 

0.10305 

0.17500 

0.24329 

0.30000 

0.32561 

lls4p 

46 

2162 

0.09343 

0.16698 

0.23127 

0.28631 

0.31915 

12s4p 

48 

2352 

0.08631 

0.15221 

0.21088 

0.26531 

0.29847 

13s4p 

50 

2550 

0.07882 

0.13922 

0.19333 

0.24588 

0.27490 

14s4p 

52 

2756 

0.07438 

0.13171 

0.18541 

0.23403 

0.27104 

14s5p 

58 

3422 

0.06517 

0.11894 

0.16774 

0.21566 

0.24839 

15s5p 

60 

3660 

0.06421 

0.11339 

0.15984 

0.20628 

0.24809 

16s5p 

62 

3906 

0 . 06042 

0.10701 

0.15182 

0.19636 

0.23374 

17s5p 

64 

4160 

0.05433 

0.10096 

0.14663 

0.18510 

0.22524 

17s6p 

70 

4970 

0.05131 

0.09598 

0.13461 

0.17445 

0.21288 

18s6p 

72 

5256 

0.04909 

0.09170 

0.13033 

0.16591 

0.20358 

"includes  orbital  products 


and  derivatives  {-r=^  <t>  } 
dRij 
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Table  5 


Least  squares  optimal  values  of  a,  and  -y  for  the 
function  r  -  atanh(/SM  )  at  cutoff  tolerances  of  10  , 

10  6 ,  10  9 ,  10  12 ,  and  1015  for  C  two -electron  integrals. 


tolerance 

a 

0 

7 

M 

O 

i 

o 

Ol 

186 

0.011401 

1.063459 

io*06 

400 

0.007137 

1.173166 

io*09 

528 

0.003003 

1.454537 

io*12 

726 

0.001181 

1.679749 

•o 

H 

< 

O 

tH 

1012 

0.000665 

1.761226 
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Table  6 


Lease  squares  optimal  values  of  a,  f),  and  7  for  the 
function  r  -  atanh(^M^)  at  cutoff  tolerances  of  10~3, 

10  8 ,  10  9 ,  10  12 ,  and  1015  for  C  two -electron  Integrals 
and  derivatives 


tolerance 

a 

0 

1 

io*03 

326 

0.013045 

1.027928 

<o 

0 

1 

0 

H 

605 

0.007877 

1 . 143445 

10  09 

772 

0.003882 

1.360329 

M 

H 

1 

O 

977 

0.001807 

1.551929 

10 


-15 


1324 


0.001100 


1.616605 
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Table  7  Fraction  of  rank  at  cutoff  tolerances  of  10  ,  10  ,  10  , 

-12  -15 

10  ,  and  10  for  two-electron  integrals  and 

integrals  plus  derivatives  for  the  Dunning  and  6-311g 
basis  sets. 


frac . of 

frac . of 

frac . of 

frac. of 

frac. of 

basis  set  # 

basis 

full 

rank  at 

rank  at 

rank  at 

rank  at 

rank  at 

function 

rank 

10**- 3 

10**- 6 

10** -9 

10**- 12 

10**- 15 

Dunning 

integrals 

6-311g 

28 

406 

0.15517 

0.34236 

0.50985 

0.56897 

0.58128 

integrals 

Dunning 

26 

351 

0.19088 

0.36752 

0.45014 

0.47009 

0.48718 

integrals  + 
derivatives 

28 

812 

0.14039 

0.25123 

0.33867 

0.36823 

0.37438 

6-311g 
integrals  + 
derivatives 

26 

702 

0.16524 

0.25783 

0.31624 

0.31197 

0.32051 
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Table  8  Fraction  of  rank  at  cutoff  tolerances  of  10  3 ,  10  6 ,  10 

-12  “15 

10  ,  and  10  for  benzene  two -electron  integrals  and 

integrals  plus  derivatives  with  an  ST0-3G  basis. 


type  of 
calc. 

#  basis 
functions 

full 

rank 

frac.of 
rank  at 
10**- 3 

frac . of 
rank  at 
10**- 6 

frac . of 
rank  at 
10** -9 

frac.of 
rank  at 
10**- 12 

frac.of 
rank  at 
10**- 15 

integrals 

36 

666 

0 . 14414 

0.31081 

0.46547 

0.58859 

0.63814 

integrals  + 
derivatives 

36 

1332 

0.15165 

0.31081 

0.46547 

0.58859 

0.63814 

Number  of  Orbital  Basis  Functions 


Figure  1.  Rank  reduction  for  diatomic  carbon  two-electron 

“3 

integrals  with  even-tempered  basis  sets  at  tolerances  of  10  ' 
10'V  10'9.  10‘12.  and  10‘13. 


Number  of  Orbital  Basis  Functions 


Figure  2.  Rank  reduction  for  diatomic  carbon  two-electron 
integrals  and  derivatives  with  even-tempered  basis  sets  at 
tolerances  of  10  3,  10  8,  10  9,  10  12 ,  and  10  15 . 


Number  of  Orbitai  Basis  Functions 


Figure  3.  Fractional  rank  r/M  for  diatomic  carbon  two -electron 
integrals  with  even-tempered  basis  sets  at  tolerances  of  10" 3 , 
10'6,  10"9,  10"12,  and  10'15. 


Fractional  Rank  of  Matrix 
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Figure  4.  Fractional  rank  r/M  for  diatomic  carbon  two-electron 
integrals  and  derivatives  with  even-tempered  basis  sets  at 
tolerances  of  10  3 ,  10  6 ,  10  9 ,  10  12 ,  and  10  13 . 


0  400  800  1200  1600  2000  2400  2800 

M  (number  of  basis  functions  squared) 


Figure  5.  Least  squares  fit  to  actual  data  for  diatomic  carbon 
two-electron  integrals  with  even-tempered  basis  sets  at  tolerances 


of  10 


-3 


10 


-8 


10 


10 


-12 


and  10 


-13 


r 

Q  L—J I I I 1 I I  1 I'll 1  1  '  1  I  ''  I  I  I  I 

0  1000  2000  3000  4000  5000  60( 

M  (number  of  basis  functions  squared) 


Figure  6.  Leasts  squares  fit  to  actual  data  for  diatomic  carbon 
two-electron  Integrals  and  derivatives  for  even-tempered  basis 
sets  at  tolerances  of  10*3,  10* 8 ,  10"9,  10"1Z,  and  10*1S. 
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Figure  7.  Ratio  of  CPU  time  weighted  arithemetic  operactions*  for 
Cholesky  algorithm  to  conventional  methods  for  diatomic  carbon 
two-electron  integrals  with  even-tempered  basis  sets  at  tolerances 
10'3,  10'8,  10*9f  10'12,  and  10'15. 


Integral  evaluation  time  (503  /jsec)  and  floating  point 
multiplication  (0.540  nsec)  are  included  in  these  timings 
(see  text). 
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Figure  8.  Redo  of  CPU  time  weighted  arithemetic  operactions*  for 
Cholesky  algorithm  to  conventional  methods  for  diatomic  carbon 
two-electron  integrals  and  derivatives  with  even-tempered  basis 
sets  at  tolerances  of  10  3 ,  10'9,  10" 9 ,  10  12 ,  and  10  15 . 


